Librerías
Las librerías utilizadas son las siguientes:
library(olsrr)
library(BSDA)
library(MVN)
library(tidyverse)
library(hrbrthemes)
library(viridis)
library(ggplot2)
library(corrplot)
library(psych)
library(lmtest)
library(car)
library(tidymodels)
library(reshape2)
library(MLmetrics)
library(plotly)Carga de Datos
data(trees)
attach(trees)
head(trees)
#> Girth Height Volume
#> 1 8.3 70 10.3
#> 2 8.6 65 10.3
#> 3 8.8 63 10.2
#> 4 10.5 72 16.4
#> 5 10.7 81 18.8
#> 6 10.8 83 19.7Preguntas
Pregunta 1
Construya gráficas apropiadas para las variables (histogramas, diagramas de dispersión por pares, boxplots para datos simétricos y asimétricos
par(mfrow=c(2,2))
boxplot(trees$Height, main='Height', col = 'red')
boxplot(trees$Girth, main='Girth', col = 'blue')
hist(trees$Height, main='Histogram of Height')
hist(trees$Girth, main='Histogram of Girth')
par(mfrow=c(1,1))
ggplot(trees, aes(x=Girth, y=Volume)) + geom_point(size=3) + geom_smooth(method = lm)
ggplot(trees, aes(x=Height, y=Volume)) + geom_point(size=3) + geom_smooth(method = lm)
ggplot(trees, aes(x=Height, y=Girth)) + geom_point(size=3)
cor(trees)
#> Girth Height Volume
#> Girth 1.0000000 0.5192801 0.9671194
#> Height 0.5192801 1.0000000 0.5982497
#> Volume 0.9671194 0.5982497 1.0000000
pairs.panels(trees,
smooth = TRUE, # Si TRUE, dibuja ajuste suavizados de tipo loess
scale = FALSE, # Si TRUE, escala la fuente al grado de correlación
density = TRUE, # Si TRUE, añade histogramas y curvas de densidad
ellipses = TRUE, # Si TRUE, dibuja elipses
method = "pearson", # Método de correlación (también "spearman" o "kendall")
pch = 21, # Símbolo pch
lm = FALSE, # Si TRUE, dibuja un ajuste lineal en lugar de un ajuste LOESS
cor = TRUE, # Si TRUE, agrega correlaciones
jiggle = FALSE, # Si TRUE, se añade ruido a los datos
factor = 2, # Nivel de ruido añadido a los datos
hist.col = 4, # Color de los histogramas
stars = TRUE, # Si TRUE, agrega el nivel de significación con estrellas
ci = TRUE) 
Pregunta 2
Ejecute la prueba de bondad de ajuste a una normal trivariada. Interprete.
mvn(trees,mvnTest = c("mardia")) # NO
#> $multivariateNormality
#> Test Statistic p value Result
#> 1 Mardia Skewness 20.9802564628517 0.0212316600572144 NO
#> 2 Mardia Kurtosis -0.163815893754234 0.869876079301746 YES
#> 3 MVN <NA> <NA> NO
#>
#> $univariateNormality
#> Test Variable Statistic p value Normality
#> 1 Shapiro-Wilk Girth 0.9412 0.0889 YES
#> 2 Shapiro-Wilk Height 0.9655 0.4034 YES
#> 3 Shapiro-Wilk Volume 0.8876 0.0036 NO
#>
#> $Descriptives
#> n Mean Std.Dev Median Min Max 25th 75th Skew Kurtosis
#> Girth 31 13.24839 3.138139 12.9 8.3 20.6 11.05 15.25 0.5010559 -0.7109412
#> Height 31 76.00000 6.371813 76.0 63.0 87.0 72.00 80.00 -0.3568773 -0.7233677
#> Volume 31 30.17097 16.437846 24.2 10.2 77.0 19.40 37.30 1.0132739 0.2460393mvn(trees,mvnTest = c("hz")) # NO
#> $multivariateNormality
#> Test HZ p value MVN
#> 1 Henze-Zirkler 0.92118 0.03216314 NO
#>
#> $univariateNormality
#> Test Variable Statistic p value Normality
#> 1 Shapiro-Wilk Girth 0.9412 0.0889 YES
#> 2 Shapiro-Wilk Height 0.9655 0.4034 YES
#> 3 Shapiro-Wilk Volume 0.8876 0.0036 NO
#>
#> $Descriptives
#> n Mean Std.Dev Median Min Max 25th 75th Skew Kurtosis
#> Girth 31 13.24839 3.138139 12.9 8.3 20.6 11.05 15.25 0.5010559 -0.7109412
#> Height 31 76.00000 6.371813 76.0 63.0 87.0 72.00 80.00 -0.3568773 -0.7233677
#> Volume 31 30.17097 16.437846 24.2 10.2 77.0 19.40 37.30 1.0132739 0.2460393mvn(trees,mvnTest = c("royston")) # NO
#> $multivariateNormality
#> Test H p value MVN
#> 1 Royston 8.331369 0.0170622 NO
#>
#> $univariateNormality
#> Test Variable Statistic p value Normality
#> 1 Shapiro-Wilk Girth 0.9412 0.0889 YES
#> 2 Shapiro-Wilk Height 0.9655 0.4034 YES
#> 3 Shapiro-Wilk Volume 0.8876 0.0036 NO
#>
#> $Descriptives
#> n Mean Std.Dev Median Min Max 25th 75th Skew Kurtosis
#> Girth 31 13.24839 3.138139 12.9 8.3 20.6 11.05 15.25 0.5010559 -0.7109412
#> Height 31 76.00000 6.371813 76.0 63.0 87.0 72.00 80.00 -0.3568773 -0.7233677
#> Volume 31 30.17097 16.437846 24.2 10.2 77.0 19.40 37.30 1.0132739 0.2460393mvn(trees,mvnTest = c("dh")) # NO
#> $multivariateNormality
#> Test E df p value MVN
#> 1 Doornik-Hansen 124.2377 6 2.096607e-24 NO
#>
#> $univariateNormality
#> Test Variable Statistic p value Normality
#> 1 Shapiro-Wilk Girth 0.9412 0.0889 YES
#> 2 Shapiro-Wilk Height 0.9655 0.4034 YES
#> 3 Shapiro-Wilk Volume 0.8876 0.0036 NO
#>
#> $Descriptives
#> n Mean Std.Dev Median Min Max 25th 75th Skew Kurtosis
#> Girth 31 13.24839 3.138139 12.9 8.3 20.6 11.05 15.25 0.5010559 -0.7109412
#> Height 31 76.00000 6.371813 76.0 63.0 87.0 72.00 80.00 -0.3568773 -0.7233677
#> Volume 31 30.17097 16.437846 24.2 10.2 77.0 19.40 37.30 1.0132739 0.2460393mvn(trees,mvnTest = c("energy")) # NO
#> $multivariateNormality
#> Test Statistic p value MVN
#> 1 E-statistic 1.191647 0.005 NO
#>
#> $univariateNormality
#> Test Variable Statistic p value Normality
#> 1 Shapiro-Wilk Girth 0.9412 0.0889 YES
#> 2 Shapiro-Wilk Height 0.9655 0.4034 YES
#> 3 Shapiro-Wilk Volume 0.8876 0.0036 NO
#>
#> $Descriptives
#> n Mean Std.Dev Median Min Max 25th 75th Skew Kurtosis
#> Girth 31 13.24839 3.138139 12.9 8.3 20.6 11.05 15.25 0.5010559 -0.7109412
#> Height 31 76.00000 6.371813 76.0 63.0 87.0 72.00 80.00 -0.3568773 -0.7233677
#> Volume 31 30.17097 16.437846 24.2 10.2 77.0 19.40 37.30 1.0132739 0.2460393Como se puede observar, ninguna de las pruebas de bondad de ajuste mencionan que es los datos son Normales Multivariados, ya sea por Mardia, Henze-Zirkler, Royston, Doornik-Hansen y E-statistic. Sin embargo, por fines académicos se seguira trabajando con una regresión lineal múltiple.
Pregunta 3
Calcule las correlaciones con la función mixedcor de R
cor(trees)
#> Girth Height Volume
#> Girth 1.0000000 0.5192801 0.9671194
#> Height 0.5192801 1.0000000 0.5982497
#> Volume 0.9671194 0.5982497 1.0000000
mixedCor(trees)
#> Call: mixedCor(data = trees)
#> Girth Heght Volum
#> Girth 1.00
#> Height 0.52 1.00
#> Volume 0.97 0.60 1.00Pregunta 4
Encuentre los estimadores puntuales y la ecuación de ajuste. Considere la variable “Volume” como la variable respuesta.
modelo <- lm(Volume~Girth + Height, data = trees)
summary(modelo)
#>
#> Call:
#> lm(formula = Volume ~ Girth + Height, data = trees)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -6.4065 -2.6493 -0.2876 2.2003 8.4847
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -57.9877 8.6382 -6.713 2.75e-07 ***
#> Girth 4.7082 0.2643 17.816 < 2e-16 ***
#> Height 0.3393 0.1302 2.607 0.0145 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 3.882 on 28 degrees of freedom
#> Multiple R-squared: 0.948, Adjusted R-squared: 0.9442
#> F-statistic: 255 on 2 and 28 DF, p-value: < 2.2e-16\(Volume = -57.9877 + 4.7082*Girth + 0.3393*Height\)
Pregunta 5
Interprete el coeficiente de regresión asociada a la variable Girth
Por cada pulgada incrementada en la variable Girth, se reflejará como un aumento de 4.7082 pies cubicos en la variable respuesta “Volume”. La variable es significativa ya que tiene un valor p menor al valor de significancia de 0.05* (2e-16)
Pregunta 6
Presente la forma de calcular el Se. Interprete.
predicts <- cbind(predict(modelo, trees[,-3]))
originals <- c(trees$Volume)
se_table <- cbind(predicts, originals)
colnames(se_table) <- c('predicts', 'originals')
se_table_melt <- melt(se_table, id.vars='index', variable.name = 'series')
head(se_table)
#> predicts originals
#> 1 4.837660 10.3
#> 2 4.553852 10.3
#> 3 4.816981 10.2
#> 4 15.874115 16.4
#> 5 19.869008 18.8
#> 6 21.018327 19.7ggplot(se_table_melt, aes(Var1, value)) + geom_line(aes(colour = Var2))
Calculo del SE:
se_table <- data.frame(se_table)
se_table$SS <- ((se_table$originals-se_table$predicts)^2)
head(se_table)
#> predicts originals SS
#> 1 4.837660 10.3 29.8371621
#> 2 4.553852 10.3 33.0182211
#> 3 4.816981 10.2 28.9768907
#> 4 15.874115 16.4 0.2765548
#> 5 19.869008 18.8 1.1427790
#> 6 21.018327 19.7 1.7379860
SE <- sqrt(sum(se_table$SS)/(length(se_table$originals)-3))
print(SE)
#> [1] 3.881832SE <- sqrt(sum(se_table$SS)/(length(se_table$originals)-3))
print(SE)
#> [1] 3.881832Este valor también se puede obtener de la siguiente manera:
summary(modelo)
#>
#> Call:
#> lm(formula = Volume ~ Girth + Height, data = trees)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -6.4065 -2.6493 -0.2876 2.2003 8.4847
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -57.9877 8.6382 -6.713 2.75e-07 ***
#> Girth 4.7082 0.2643 17.816 < 2e-16 ***
#> Height 0.3393 0.1302 2.607 0.0145 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 3.882 on 28 degrees of freedom
#> Multiple R-squared: 0.948, Adjusted R-squared: 0.9442
#> F-statistic: 255 on 2 and 28 DF, p-value: < 2.2e-16Residual standard error : \(3.882 = SE\)
Pregunta 7
Construya intervalos de confianza para el promedio del volumen y para un valor individual para Girth=11.5 y Height=73.
y <- mean(trees$Volume)
# promedio
tsum.test(mean.x=y, s.x=SE, n.x=length(se_table$originals), conf.level=0.95, alternative = 'two.sided')$conf
#> [1] 28.74710 31.59484
#> attr(,"conf.level")
#> [1] 0.95
# valor en especifico
new <- data.frame(cbind(11.5, 73))
predict(modelo, newdata = new, interval = 'prediction', level = 0.95)
#> fit lwr upr
#> 1 4.837660 -3.561809 13.23713
#> 2 4.553852 -3.962908 13.07061
#> 3 4.816981 -3.809144 13.44311
#> 4 15.874115 7.690594 24.05764
#> 5 19.869008 11.451358 28.28666
#> 6 21.018327 12.469920 29.56673
#> 7 16.192688 7.797083 24.58829
#> 8 19.245949 11.092268 27.39963
#> 9 21.413021 13.103698 29.72234
#> 10 20.187581 12.047515 28.32765
#> 11 22.015402 13.775543 30.25526
#> 12 21.468465 13.327934 29.60900
#> 13 21.468465 13.327934 29.60900
#> 14 20.506154 12.270386 28.74192
#> 15 23.954110 15.854249 32.05397
#> 16 27.852203 19.760075 35.94433
#> 17 31.583966 23.126434 40.04150
#> 18 33.806482 25.303646 42.30932
#> 19 30.600978 22.388655 38.81330
#> 20 28.697035 19.945837 37.44823
#> 21 34.388184 26.295494 42.48087
#> 22 36.008319 27.878179 44.13846
#> 23 35.385260 27.237521 43.53300
#> 24 41.768998 33.386120 50.15188
#> 25 44.877702 36.655198 53.10021
#> 26 50.942868 42.647210 59.23853
#> 27 52.223751 43.899133 60.54837
#> 28 53.428513 45.064545 61.79248
#> 29 53.899329 45.522483 62.27618
#> 30 53.899329 45.522483 62.27618
#> 31 68.515305 59.707135 77.32347Pregunta 8
Es la variable Height significativa? Ejecute la prueba de hipótesis completa.Presente la forma de calcular el p-valor.
H0 = No es significativo H1 = Es significtivo
anova(modelo)
#> Analysis of Variance Table
#>
#> Response: Volume
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Girth 1 7581.8 7581.8 503.1503 < 2e-16 ***
#> Height 1 102.4 102.4 6.7943 0.01449 *
#> Residuals 28 421.9 15.1
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
2*pt(2.607, 28, lower.tail = FALSE)
#> [1] 0.01447721Dado los resultados, se puede decir que la variable Height es significativa, ya que el valor p es menor al valor de alpha (0.05) por esto no se rechaza la hipótesis nula
Pregunta 9
Es el modelo significativo? Ejecute la prueba de hipótesis completa.Presente la forma de calcular el p-valor.
H0 = No es significativo H1 = Es significtivo
g <- summary(modelo)
g$fstatistic
#> value numdf dendf
#> 254.9723 2.0000 28.0000
pf(254.9723,3,28, lower.tail = FALSE)
#> [1] 1.998157e-20El modelo es significativo ya que el valor p calculado es menor al alpha establecido, por ende no se rechaza la hipótesis nula.
Pregunta 10
Presente la forma de calcular el R2 ajustado e interprete
El valor de R2 se puede obtener mediante el siguiente comando:
summary(modelo)
#>
#> Call:
#> lm(formula = Volume ~ Girth + Height, data = trees)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -6.4065 -2.6493 -0.2876 2.2003 8.4847
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -57.9877 8.6382 -6.713 2.75e-07 ***
#> Girth 4.7082 0.2643 17.816 < 2e-16 ***
#> Height 0.3393 0.1302 2.607 0.0145 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 3.882 on 28 degrees of freedom
#> Multiple R-squared: 0.948, Adjusted R-squared: 0.9442
#> F-statistic: 255 on 2 and 28 DF, p-value: < 2.2e-16Pero también se puede calcular de la siguiente forma
r2 <- anova(modelo)
r2 <-1 - (r2$`Sum Sq`[3]/(sum(r2$`Sum Sq`)))
r2_adj <- 1 - (((1-r2)*(length(trees$Volume)-1))/(length(trees$Volume)-length(trees[,-3])-1))
print(r2)
#> [1] 0.94795
print(r2_adj)
#> [1] 0.9442322Como resultado se tiene que el modelo tiene un 94.8% de que los datos se ajusten al modelo establecido con las 2 variables (Height y Girth) y así obtener el Volume
trees$predicted <- predict(modelo)
trees$residuals <- residuals(modelo)
ggplot(trees, aes(x = Height, y = Volume)) +
geom_smooth(method = "lm", se = FALSE, color = "lightgrey") + # regression line
geom_segment(aes(xend = Height, yend = predicts), alpha = .2) + # draw line from point to line
geom_point(aes(color = abs(residuals), size = abs(residuals))) + # size of the points
scale_color_continuous(low = "green", high = "red") + # colour of the points mapped to residual size - green smaller, red larger
guides(color = FALSE, size = FALSE) + # Size legend removed
geom_point(aes(y = predicts), shape = 1) +
theme_bw()
ggplot(trees, aes(x = Girth, y = Volume)) +
geom_smooth(method = "lm", se = FALSE, color = "lightgrey") + # regression line
geom_segment(aes(xend = Girth, yend = predicts), alpha = .2) + # draw line from point to line
geom_point(aes(color = abs(residuals), size = abs(residuals))) + # size of the points
scale_color_continuous(low = "green", high = "red") + # colour of the points mapped to residual size - green smaller, red larger
guides(color = FALSE, size = FALSE) + # Size legend removed
geom_point(aes(y = predicts), shape = 1) +
theme_bw()
Pregunta 11
Calcule los residuales. Interprete.
predicts <- cbind(predict(modelo, trees[,-3]))
originals <- c(trees$Volume)
se_table <- cbind(predicts, originals)
colnames(se_table) <- c('predicts', 'originals')
se_table <- data.frame(se_table)
se_table$residuals <- residuals(modelo)
head(se_table)
#> predicts originals residuals
#> 1 4.837660 10.3 5.4623403
#> 2 4.553852 10.3 5.7461484
#> 3 4.816981 10.2 5.3830187
#> 4 15.874115 16.4 0.5258848
#> 5 19.869008 18.8 -1.0690084
#> 6 21.018327 19.7 -1.3183270El como se puede observar, existe un grado de cuanto a la predicción del Volume con el modelo, principalmente con los 3 primeros datos, pudiendo ser estos outliers. Para obtener un mejor modelo se recomienda poder extraer una mayor cantidad de datos.
H0 = Los residuales son homocedasticos, varianza constante H1 = Los residuales no homocedasticos (Heterosedasticos)
Graficar los residuos frente a los valores ajustados por el modelo, e identificar si existe un patrón cónico u otro patrón. Idealmente deberían distribuirse de forma aleatoria en torno a 0. También podemos recurrir al test de Breusch-Pagan como contraste de homocedasaticidad.
#no estandarizada
plot(modelo, which=1, col=c("red"))
bptest(modelo)
#>
#> studentized Breusch-Pagan test
#>
#> data: modelo
#> BP = 2.4681, df = 2, p-value = 0.2911
valor p = 0.2911 > 0.05, no se rechaza la H0, es decir que los residuos tienen varianza constante
Pregunta 12
Calcule los residuales estandarizados y estudentizados. Interprete.
#estandarizada y estudientizado
se_table$residuals_standard <- rstandard(modelo)
se_table$residuals_studentized <- rstudent(modelo)
# standarized graph
plot(modelo, which=2, col=c("red"))
ols_plot_resid_stand(modelo)
# studentized graph
hist(se_table$residuals_studentized, freq = FALSE)
qqPlot(modelo)
#> [1] 18 31
ols_plot_resid_stud(modelo)
Lo que puede observar es que los valores residuales estandarizados y estudiantizados tienden a acercarse a la linea de 45º, por lo que, podemos decir que los valores de los residuales poseen una distribución normal
Pregunta 13
Presente la forma de calcular la distancia de Cook de los cinco primeros registros. Interprete.
ols_plot_cooksd_bar(modelo)
ols_plot_resid_lev(modelo)
Se puede observar que los registros 20 y 31 son outliers
cooksD <- cooks.distance(modelo)
se_table$cooksD <- cbind(cooksD)
head(se_table)
#> predicts originals residuals residuals_standard residuals_studentized
#> 1 4.837660 10.3 5.4623403 1.4964901 1.5320694
#> 2 4.553852 10.3 5.7461484 1.6029462 1.6516683
#> 3 4.816981 10.2 5.3830187 1.5284555 1.5677398
#> 4 15.874115 16.4 0.5258848 0.1396700 0.1372010
#> 5 19.869008 18.8 -1.0690084 -0.2936751 -0.2888284
#> 6 21.018327 19.7 -1.3183270 -0.3696163 -0.3638447
#> cooksD
#> 1 0.0977927647
#> 2 0.1478462779
#> 3 0.1673192079
#> 4 0.0004091116
#> 5 0.0039449244
#> 6 0.0084012068influential <- cooksD[(cooksD > (3*mean(cooksD, na.rm = TRUE)))]
influential
#> 3 18 31
#> 0.1673192 0.1775359 0.6052326cooks_table <- se_tableCalculo manual:
Hmat <- hatvalues(modelo)
cooks_table$h <- c(Hmat)
data(trees)
attach(trees)
head(trees)
#> Girth Height Volume
#> 1 8.3 70 10.3
#> 2 8.6 65 10.3
#> 3 8.8 63 10.2
#> 4 10.5 72 16.4
#> 5 10.7 81 18.8
#> 6 10.8 83 19.7
k <- length(trees[,-3])
cooks_table$cookD_calc <- ((cooks_table$residuals_standard^2)*cooks_table$h)/((k+1)*(1-cooks_table$h)) # k+1 => k: numero de variables predictorasEn la tabla cooks_table, en la columna “cookD_calc” se muestran los resultados de la distancia de cook que antes se calculo con la formula:
head(cooks_table, 5)
#> predicts originals residuals residuals_standard residuals_studentized
#> 1 4.837660 10.3 5.4623403 1.4964901 1.5320694
#> 2 4.553852 10.3 5.7461484 1.6029462 1.6516683
#> 3 4.816981 10.2 5.3830187 1.5284555 1.5677398
#> 4 15.874115 16.4 0.5258848 0.1396700 0.1372010
#> 5 19.869008 18.8 -1.0690084 -0.2936751 -0.2888284
#> cooksD h cookD_calc
#> 1 0.0977927647 0.11582883 0.0977927647
#> 2 0.1478462779 0.14720958 0.1478462779
#> 3 0.1673192079 0.17686186 0.1673192079
#> 4 0.0004091116 0.05919131 0.0004091116
#> 5 0.0039449244 0.12066468 0.0039449244Presente detalladamente la forma de calcular la matriz variancias-covariancias del vector de estimadores de beta.
matrix_trees <- trees[, 1:2]
matrix_trees$r <- c(rep(1, 31))
matrix_trees <- matrix_trees[, c(3,1,2)]
matrix_trees <- as.matrix(matrix_trees)
vcov(modelo)
#> (Intercept) Girth Height
#> (Intercept) 74.6189461 0.43217138 -1.05076889
#> Girth 0.4321714 0.06983578 -0.01786030
#> Height -1.0507689 -0.01786030 0.01693933
(SE^2)*solve(t(matrix_trees)%*%(matrix_trees))
#> r Girth Height
#> r 74.6189461 0.43217138 -1.05076889
#> Girth 0.4321714 0.06983578 -0.01786030
#> Height -1.0507689 -0.01786030 0.01693933Pregunta 15
Calcule el VIF. Comente sus resultados
vif(modelo)
#> Girth Height
#> 1.36921 1.36921El factor de inflación de varianza (VIF) muestra que tanto para la variable Girth como Height toma valores menores a 10, lo que significa que ambas variables entran al modelo.
Pregunta 16
Elija al azar 70% de los datos y proceda a ejecutar la RLM para obtener la ecuación de ajuste.
division <- trees %>% initial_split(prop = 0.7, strata = Volume)
training_set <- division %>% training()
testing_set <- division %>% testing()
length(training_set$Volume)
#> [1] 21
length(testing_set$Volume)
#> [1] 10Desarrollo de RLM:
modelo_train <- lm(Volume ~ Girth + Height, data = training_set)
summary(modelo_train)
#>
#> Call:
#> lm(formula = Volume ~ Girth + Height, data = training_set)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -6.1422 -2.1230 -0.3601 2.9643 5.1460
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -48.9374 13.2267 -3.700 0.00164 **
#> Girth 4.6903 0.3118 15.043 1.23e-11 ***
#> Height 0.2267 0.1885 1.203 0.24469
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 3.693 on 18 degrees of freedom
#> Multiple R-squared: 0.9447, Adjusted R-squared: 0.9385
#> F-statistic: 153.6 on 2 and 18 DF, p-value: 4.864e-12Pregunta 17
Con la ecuación obtenida en el paso anterior, calcule los valores ajustados del 30% restante
predicts <- cbind(predict(modelo, testing_set[, 1:2]))
predicts_table <- cbind(predicts, c(testing_set[,3]))
colnames(predicts_table) <- c('predicts', 'originals')
predicts_table <- as.data.frame(predicts_table)
predicts_table
#> predicts originals
#> 3 4.816981 10.2
#> 4 15.874115 16.4
#> 6 21.018327 19.7
#> 7 16.192688 15.6
#> 14 20.506154 21.3
#> 19 30.600978 25.7
#> 20 28.697035 24.9
#> 22 36.008319 31.7
#> 30 53.899329 51.0
#> 31 68.515305 77.0Pregunta 18
Calcule el MSE del 30% restante.
mean((predicts_table$original-predicts_table$predicts)^2)
#> [1] 16.93677
MSE(predicts_table$predicts, predicts_table$original)
#> [1] 16.93677Pregunta 19
Explique la diferencia entre el cuadrado medio del error del análisis de variancia y el MSE
Según los resultados del anova el MSE es 12.1 (es el resultado de dividir la suma de cuadrados del error entre n - k - 1) donde (k es numero de coeficientes = 3, n es el tamaño de la muestra = 21), mientras que el MSE calculado resultó 13.47 (es el resultado de la división entre n).
Pregunta 20
Grafique los datos en 3D usando R
fig <- plot_ly(x = trees$Girth, y = trees$Height, z = trees$Volume)
fig